#include "sgPardiso.h"
#include "gtest/gtest.h"

//#include "DataStructure/Common/include/Constant.h"

#include <complex.h>

using namespace SG::Algebra;  // Open the namespace
//using SG::DataStructure::Common::TOLERANCE;

class PardisoTest : public ::testing::Test
{
  public:
    static void SetUpTestSuite ()
    {   
    }
};

TEST_F (PardisoTest, singleTest)
{
    TPardiso<Real_t> solver;

    SpMat_t mat;

    typedef Eigen::Triplet<Real_t> Triplet;
    std::vector<Triplet>           triplets;

    int m = 1;  // row count
    int n = 1;  // column count

    triplets.push_back (Triplet (0, 0, 6.0));

    // Algebra::SpMat_t mat (m, n);
    mat.resize (m, n);
    mat.setFromTriplets (triplets.begin (), triplets.end ());

    solver.compute (mat);

    // Check
    auto info1 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info1);

    // solution from OctaveOnline
    // x = inv(A)*b
    // see https://octave-online.net/
    Real_t exact[] = { 1 };

    Vec_t b (m);
    b[0] = 6.0;

    Vec_t x (m);
    solver.solve (b, x);

    // Check
    auto info2 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info2);

    // for (int i = 0; i < m; ++i)
    // {
    //     std::cout << exact[i] << ", " << x[i] << std::endl;
    // }

    for (int i = 0; i < m; ++i)
    {
        ASSERT_NEAR (exact[i], x (i), TOLERANCE);
    }
}

/*
 * Test the A*x = b given by SuperLU User Guide.
 * where A is a 3X3 matrix and b is a 3X1 matrix.
 *
 * A =
 *   [
 *       06.00,  07.00,  05.00;
 *       07.00,  13.00,  08.00;
 *       05.00,  08.00,  06.00
 *   ]
 *
 * b = [9.0; 10.0; 9.0]
 *
 * OctaveOnline give the solution as followings,
 * x = inv(A)*b = [1; -1; 2]
 *
 * see https://octave-online.net/
 */
TEST_F (PardisoTest, realTest)
{
    TPardiso<Real_t> solver;

    SpMat_t mat;

    typedef Eigen::Triplet<Real_t> Triplet;
    std::vector<Triplet>           triplets;

    int m = 3;  // row count
    int n = 3;  // column count

    triplets.push_back (Triplet (0, 0, 6.0));
    triplets.push_back (Triplet (0, 1, 7.0));
    triplets.push_back (Triplet (0, 2, 5.0));

    triplets.push_back (Triplet (1, 0, 7.0));
    triplets.push_back (Triplet (1, 1, 13.0));
    triplets.push_back (Triplet (1, 2, 8.0));

    triplets.push_back (Triplet (2, 0, 5.0));
    triplets.push_back (Triplet (2, 1, 8.0));
    triplets.push_back (Triplet (2, 2, 6.0));

    // Algebra::SpMat_t mat (m, n);
    mat.resize (m, n);
    mat.setFromTriplets (triplets.begin (), triplets.end ());

    solver.compute (mat);

    // Check
    auto info1 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info1);

    // solution from OctaveOnline
    // x = inv(A)*b
    // see https://octave-online.net/
    Real_t exact[] = { 1, -1, 2 };

    Vec_t b (m);
    b[0] = 9.0;
    b[1] = 10.0;
    b[2] = 9.0;

    Vec_t x (m);
    solver.solve (b, x);

    // Check
    auto info2 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info2);

    // for (int i = 0; i < m; ++i)
    // {
    //     std::cout << exact[i] << ", " << x[i] << std::endl;
    // }

    for (int i = 0; i < m; ++i)
    {
        ASSERT_NEAR (exact[i], x (i), TOLERANCE);
    }
}

TEST_F (PardisoTest, complexTest)
{
    TPardiso<Complex_t> solver;

    ComplexSpMat_t mat;

    typedef Eigen::Triplet<Complex_t> Triplet;
    std::vector<Triplet>           triplets;

    int m = 3;  // row count
    int n = 3;  // column count

    triplets.push_back (Triplet (0, 0, 6.0I));
    triplets.push_back (Triplet (0, 1, 7.0I));
    triplets.push_back (Triplet (0, 2, 5.0I));

    triplets.push_back (Triplet (1, 0, 7.0I));
    triplets.push_back (Triplet (1, 1, 13.0I));
    triplets.push_back (Triplet (1, 2, 8.0I));

    triplets.push_back (Triplet (2, 0, 5.0I));
    triplets.push_back (Triplet (2, 1, 8.0I));
    triplets.push_back (Triplet (2, 2, 6.0I));

    // Algebra::SpMat_t mat (m, n);
    mat.resize (m, n);
    mat.setFromTriplets (triplets.begin (), triplets.end ());

    solver.compute (mat);

    // Check
    auto info1 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info1);

    // solution from OctaveOnline
    // x = inv(A)*b
    // see https://octave-online.net/
    Complex_t exact[] = { 1, -1, 2 };

    ComplexVec_t b (m);
    b[0] = Complex_t(9.0I);
    b[1] = Complex_t(10.0I);
    b[2] = Complex_t(9.0I);

    ComplexVec_t x (m);
    solver.solve (b, x);

    // Check
    auto info2 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info2);

    // for (int i = 0; i < m; ++i)
    // {
    //     std::cout << exact[i] << ", " << x[i] << std::endl;
    // }

    for (int i = 0; i < m; ++i)
    {
        ASSERT_NEAR (exact[i].real(), x (i).real(), TOLERANCE);
        ASSERT_NEAR (exact[i].imag(), x (i).imag(), TOLERANCE);
    }
}

TEST_F(PardisoTest, sloveInversers) {
    TPardiso<Real_t> solver;
    const int N = 6;

    // ==== Construct an upper triangular matrix ====
    SpMat_t A(N, N);
    std::vector<Eigen::Triplet<Real_t>> triplets;
    for (int i = 0; i < N; ++i) {
        triplets.emplace_back(i, i, 2.0 + i);
        for (int j = i + 1; j < N; ++j) {
            triplets.emplace_back(i, j, -0.5 * (j - i));
        }
    }
    A.setFromTriplets(triplets.begin(), triplets.end());
    A.makeCompressed();

    // ==== Construct an identity matrix ====
    SpMat_t B(N, N);
    B.setIdentity();

    // ==== Solving the inverse matrix ====
    SpMat_t invA;
    solver.compute(A);
    solver.solve(B, invA);

    // ==== Verify the product matrix ====
    Eigen::MatrixXd product = A.toDense() * invA.toDense();

    // Verification results
    Eigen::MatrixXd identity = Eigen::MatrixXd::Identity(N, N);
    double tolerance = 1e-6;
    EXPECT_TRUE(product.isApprox(identity, tolerance))
        << "Verification failed: maximum deviation=" 
        << (product - identity).cwiseAbs().maxCoeff();
}